#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <stdbool.h>
#include <quadmath.h>
#include <float.h>
#include <immintrin.h>

void check_double_range()
{
   double x;
   int i;
   char buff[64];

   printf("\n\n");
   printf("\nrun %s()\n", __func__);

   x=0x1.0p+1020;
   i=1020;
   while (1)
   {
      sprintf(buff,"2^%d\t=%la",i,x);
      x=x*2; i++;
      if ( isinf(x))
         break;
   }
   printf("%s\n",buff);

   x=0x1.0p-1020;
   i=-1020;
   while (1)
   {
      sprintf(buff,"2^%d\t=%la\n",i,x);
      x=x*0.5; i--;
      if ( x == 0.0)
         break;
   }
   printf("%s\n\n",buff);
}

void check_exp_para_range()
{
   char numbuf[128];
   char outbuf1[64];
   char outbuf2[64];
   char outbuf3[64];

   printf("\n\n");
   printf("\nrun %s()\n", __func__);

   double f= 1024- (1.0/256.0);
   printf("1024.0-(1.0/256.0)=%.16lf=%.14la\n",f,f);

   double es[]= {-1074, -1022, -960, 1022,1023,1024-1.0/256, 1024 };

   __float128 q2 = 2.0Q;
   for (int i=0;i<sizeof(es)/sizeof(es[0]);i++ )
   {
      __float128 e = es[i];
      __float128 y = powq(q2, e);
      __float128 x = logq(y);
      quadmath_snprintf (outbuf1, sizeof outbuf1, "%.20Qf", 64,x);
      quadmath_snprintf (outbuf2, sizeof outbuf2, "%.14Qa", 64,x);
      printf("exp(%s)=exp(%s)=\t2^%.12lf\n",outbuf1,outbuf2,es[i]);
   }
   printf("\n");

   double bound[]   =   {-960,                1023 };
   double bound_hi[]=   {-0x1.4cb5ecf0bp+9,   0x1.628b76e3ap+9 };

   q2 = 2.0Q;
   for (int i=0;i<sizeof(bound)/sizeof(bound[0]);i++ )
   {
      __float128 e = bound[i];
      __float128 y = powq(q2, e);
      __float128 x = logq(y);
      __float128 x_hi = bound_hi[i];
      __float128 x_low = x - x_hi;

      quadmath_snprintf (outbuf1, sizeof outbuf1, "%.20Qf", 64,x);
      quadmath_snprintf (outbuf2, sizeof outbuf2, "%.13Qa", 64,x_hi);
      quadmath_snprintf (outbuf3, sizeof outbuf3, "%.13Qa", 64,x_low);
      printf("exp(%s)=\texp(%s + %s)=\t2^%.12lf\n",outbuf1,outbuf2,outbuf3,bound[i]);
   }
   printf("\n");
}

void print_double_raw_data(double x)
{
   unsigned char *p=(unsigned char*)(&x);

   printf("%.13la=0x",x);
   for (int i=sizeof(x)-1;i>=0;i--)
      printf("%02x",p[i]);
   printf("\n");
}

double truncate(double x, int keep_bits)
{
   int e;
   frexp(x,&e);
   double scale=1.0;
   int    shift=0;
   if ( e-1 < keep_bits)
   {
      shift= keep_bits-(e-1);
      scale= ldexp(0.5, shift+1);
   }
   return floor(x*scale)/scale;
}

void get_const()
{
   printf("\nrun %s()\n", __func__);

   char outbuf[64];
   double f;
   __float128 ln2hi,ln2low;
   __float128 q2 = 2;
   __float128 q1 = 1;
   __float128 ln2 = logq(q2);

   quadmath_snprintf (outbuf, sizeof outbuf, "%.24Qa", 64,ln2);
   printf("ln(2)=%s\n",outbuf);

   quadmath_snprintf (outbuf, sizeof outbuf, "%.29Qf", 64,ln2);
   printf("ln(2)=%s\n",outbuf);

   ln2hi=0x1.62e42fefap-1q;   printf("       %.18f=%.13la\n",f,f);
   ln2low= ln2-ln2hi;

   quadmath_snprintf (outbuf, sizeof outbuf, "%.13Qa", 64,ln2hi);
   printf("ln2hi=%s\n",outbuf);

   quadmath_snprintf (outbuf, sizeof outbuf, "%.13Qa", 64,ln2low);
   printf("ln2low=%s\n",outbuf);

   __float128 rln2 = q1 / ln2;
   quadmath_snprintf (outbuf, sizeof outbuf, "%.14Qa", 64,rln2);
   printf("1/ln(2)=%s\n",outbuf);
   printf("1/ln(2)=%s\n","\"1.44269504088896340735\"");
   f=0x1.71547652b82fep+0;  printf("         %.18f=%.13la\n",f,f);
   f=0x1.71547652b82ffp+0;  printf("         %.18f=%.13la\n",f,f);

   //exp(0x1.6232bdd7abcd24p+9)=       2^1022.000000000000
   f=0x1.6232bdd7abcd2p+9;
   print_double_raw_data(f);

   f=-708.396418532264;
   print_double_raw_data(f);

   f= 709.7800052872;
   print_double_raw_data(f);
}

void test_floor_ceil()
{
   printf("\nrun %s()\n", __func__);
   double ary[]= { -1.51, -1.5, -1.49, -0.49, 0.5, 0.51, 1.49, 1.5, 1.51 };
   for (int i=0; i<sizeof(ary)/sizeof(ary[0]);i++)
   {
      double x=ary[i];
      printf("floor(%lf)=%lf\n",x,floor(x));
      printf("ceil(%lf)=%lf\n",x,ceil(x));
      printf("round(%lf)=%lf\n",x,round(x));
   }
}


void test_truncate()
{
   printf("\nrun %s()\n", __func__);
   double ary[]= { 1.324345, 132.4345 };

   for (int i=0; i<sizeof(ary)/sizeof(ary[0]);i++)
   {
      double x= ary[i];
      double x0= truncate(x,36);
      printf("truncate(%la,%d)=%la\n",x,36,x0);
   }
}

int main(int argc, char* argv[])
{
   check_double_range();
   check_exp_para_range();

   //test_floor_ceil();
   test_truncate();
   get_const();
   return 0;
}
/*
  双精度浮点数d的范围是: [2^-1074,2^1024), 左开右闭区间
  let exp(x) is double, then  <= x <=
*/

